//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-
/*************************************************************************
 * template<class T>
 * class BoostLapackSymSolver
 *
 * Input Matrices = Lower Triangle Packed (in fortran)
 * using triangle_matrix<T,upper> 
 *
 * Expert driver to calculate M eigenstates, M < N 
 *************************************************************************/
#ifndef OHMMS_OHMMSMATRIX_LAPACK_SYMSOLVER_PACKED_H
#define OHMMS_OHMMSMATRIX_LAPACK_SYMSOLVER_PACKED_H

#ifndef EIGENSPACEDIMENSION_H
#include "EigSolver/EigenSpaceDimension.h"
#endif

#ifdef OHMMS_MTL
#include "mtl/mtl.h"
#else
#include "OhmmsPETE/OhmmsVector.h"
#endif

#ifdef ADD_
#define dspevx dspevx_
#define zhpevx zhpevx_
#define dspgvx dspgvx_
#define zhpgvx zhpgvx_
#endif
extern "C" {

  void dspevx(const char& jobz, const char& range, const char& uplo, 
              const int& n, double* a, const double& vl, const double& vu,
              const int& il, const int& iu, const double& tol, const int& m, 
              double* w, double* z, const int& ldz, double* workd, 
              int* iwork, int* ifail, int& info); 

  void zhpevx(const char& jobz, const char& range, const char& uplo, 
              const int& n, complex<double>* a, 
              const double& vl, const double& vu,
              const int& il, const int& iu, const double& tol, const int& m, 
              double* w, complex<double>* z, const int& ldz, 
              complex<double>* workd, double* rworkd, int* iwork,
              int* ifail, int& info); 

  void dspgvx(const int& itype, const char& jobz, 
              const char& range, const char& uplo, 
              const int& n, double* a, double* b, 
              const double& vl, const double& vu,
              const int& il, const int& iu, const double& tol, const int& m, 
              double* w, double* z, const int& ldz, double* workd, 
              int* iwork, int* ifail, int& info); 

  void zhpgvx(const int& itype, const char& jobz, 
              const char& range, const char& uplo, 
              const int& n, complex<double>* a, complex<double>* b,
              const double& vl, const double& vu,
              const int& il, const int& iu, const double& tol, const int& m, 
              double* w, complex<double>* z, const int& ldz, 
              complex<double>* workd, double* rworkd, int* iwork,
              int* ifail, int& info); 
}

template<class T, class Ham_t, class EigVec_t>
class LapackPackedSymSolver {

public:

  typedef typename OHMMS::type_traits<T>::value_type value_type;
  typedef typename OHMMS::type_traits<T>::real_type  real_type;
  typedef Ham_t    HamiltonianMatrix_t;
  typedef EigVec_t EigenStates_t;
#ifdef OHMMS_MTL
  typedef mtl::dense1D<value_type> WorkArray_t;
  typedef mtl::dense1D<real_type>  RealArray_t;
  typedef mtl::dense1D<int>        IndexArray_t;
#else
  typedef Vector<value_type> WorkArray_t;
  typedef Vector<real_type>  RealArray_t;
  typedef Vector<int>        IndexArray_t;
#endif
  typedef typename WorkArray_t::size_type   size_type;
  EigenSpaceDimension& States;       
  HamiltonianMatrix_t  H, S;
  
  LapackPackedSymSolver(EigenSpaceDimension& dim, bool nonortho=false): 
  States(dim), Itype(1), Range('I'), Uplo('L'), Jobz('N'),MyTol(0){ 
#ifdef OHMMS_MTL
    // resize for packed storage is not working properly
    H = HamiltonianMatrix_t(States.N,States.N);
    if(nonortho) S = HamiltonianMatrix_t(States.N, States.N);
#else
    H.resize(States.N, States.N);
    if(nonortho) S.resize(States.N, States.N);
#endif
  }


  ~LapackPackedSymSolver() { }

  //bool getEigvectors(HamiltonianMatrix_t& H, EigenStates_t& eigvec) {
  template<class InHamiltonian_t>
  bool getEigvectors(InHamiltonian_t& h, EigenStates_t& Psi) {
    copy(h,H);  // H = h*phase
    Jobz = 'V';
    return solve(H.data(), Psi.data()); 
  }

  template<class InHamiltonian_t>
  bool getEigvalues(InHamiltonian_t& h) {
    copy(h,H);  // H = h*phase
    Jobz = 'N';
    return solve(H.data(), H.data()); 
  }


  template<class InHamiltonian_t, class IterType>
  bool 
  getEigvectors(InHamiltonian_t& h, IterType phi, EigenStates_t& Psi) {
    copy(h,phi,H);  // H = h*phase
    Jobz = 'V';
    return solve(H.data(), Psi.data()); 
  }

  template<class InHamiltonian_t, class IterType>
  bool getEigvalues(InHamiltonian_t& h, IterType phi) {
    copy(h,phi,H);  // H = h*phase
    Jobz = 'N';
    return solve(H.data(), H.data()); 
  }


  template<class InHamiltonian_t>
  bool 
  getEigvectors(InHamiltonian_t& h, InHamiltonian_t& s, EigenStates_t& Psi) {
    copy(h,H);  // H = h*phase
    copy(s,S);  // H = h*phase

    Jobz = 'V';
    return solve(H.data(), S.data(), Psi.data());
  }

  template<class InHamiltonian_t>
  bool getEigvalues(InHamiltonian_t& h, InHamiltonian_t& s) {
    copy(h,H);  // H = h*phase
    copy(s,S);  // H = h*phase

    Jobz = 'N';
    return solve(H.data(), S.data(), H.data());
  }


  template<class InHamiltonian_t, class IterType>
  bool 
  getEigvectors(const InHamiltonian_t& h, const InHamiltonian_t& s, 
		IterType phi, EigenStates_t& Psi) {
    OHMMS::copy(h,phi,H);  // H = h*phase
    OHMMS::copy(s,phi,S);  // H = h*phase

    Jobz = 'V';
    return solve(H.data(), S.data(), Psi.data()); 
  }

  template<class InHamiltonian_t, class IterType>
  bool 
  getEigvalues(InHamiltonian_t& h, InHamiltonian_t& s, IterType phi) {
    OHMMS::copy(h,phi,H);  // H = h*phase
    OHMMS::copy(s,phi,S);  // H = h*phase
    Jobz = 'N';
    return solve(H.data(), S.data(), H.data()); 
  }

//    template<class OtherMat_t>
//    bool copyEigvectors(OtherMat_t& psi0) {
//      copy(Psi,psi0);
//    }

  inline real_type eigval(size_type i) const { return eigVal[i];}

  int Itype;
  int Info;
  char Range, Jobz, Uplo;
  real_type MyTol;

  // \todo make these arrays static workspace
  IndexArray_t Ifail, Iwork;
  WorkArray_t workd; 
  RealArray_t rworkd;
  RealArray_t eigVal;

  //inline bool solve(double* a, double *z) {

  inline bool solve(double*  a, double *z) {
    check_workspace(States.N, 8, 0);
    dspevx(Jobz, Range, Uplo, States.N, a, 
           States.MinE, States.MaxE, States.First, States.Last, 
           MyTol, States.M, eigVal.data(), z, States.N, 
           workd.data(), Iwork.data(), Ifail.data(), Info);
    return true;
  }

  inline bool solve(double* ap, double *bp, double* z) {

    check_workspace(States.N, 8, 0);
    dspgvx(Itype, Jobz, Range, Uplo, States.N, ap, bp, 
           States.MinE, States.MaxE, States.First, States.Last, 
           MyTol, States.M, eigVal.data(), z, States.N, 
           workd.data(), Iwork.data(), Ifail.data(), Info);
    return true;
  }

  inline bool solve(complex<double>* a, complex<double> *z) {

    check_workspace(States.N, 2, 7);
    zhpevx(Jobz, Range, Uplo, States.N, a, 
           States.MinE, States.MaxE, States.First, States.Last, 
           MyTol, States.M, eigVal.data(), z, States.N, 
           workd.data(), rworkd.data(), Iwork.data(), Ifail.data(), Info);
    return true;
  }

  inline bool solve(complex<double>* ap, complex<double>* bp, complex<double>* z) {

    check_workspace(States.N, 2, 7);
    zhpgvx(Itype, Jobz, Range, Uplo, States.N, ap, bp, 
           States.MinE, States.MaxE, States.First, States.Last, 
           MyTol, States.M, eigVal.data(), z, States.N, 
           workd.data(), rworkd.data(), Iwork.data(), Ifail.data(), Info);
    return true;
  }

  void check_workspace(int n, int workd_size, int rworkd_size) {

    if(Ifail.size() < n)  Ifail.resize(n);
    if(Iwork.size() < 5*n) Iwork.resize(5*n);
    if(eigVal.size() < n) eigVal.resize(n);
    if(workd.size() < workd_size*n) workd.resize(workd_size*n);
    if(rworkd.size() < rworkd_size*n) rworkd.resize(rworkd_size*n);
  }
};

#endif

/***************************************************************************
 * $RCSfile: LapackSymSolverPacked.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: LapackSymSolverPacked.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

